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The mixing of materials due to the Richtmyer-Meshkov instability and the ensuing turbulent 
behavior is of intense interest in a variety of physical systems including inertial confinement fusion, 
combustion, and the final stages of stellar evolution. Extensive numerical and laboratory studies of 
shock-driven mixing have demonstrated the rich behavior associated with the onset of turbulence 
due to the shocks. Here we report on progress in understanding shock-driven mixing at interfaces 
between fluids of differing densities through 3D numerical simulations using the RAGE code in the 
implicit large eddy simulation context. We consider a shock tube configuration with a band of 
high density gas (SEe) embedded in low density gas (air). Shocks with a Mach number of 1.26 are 
passed through SEe bands, resulting in transition to turbulence driven by the Richtmyer-Meshkov 
instability. The system is followed as a rarefaction wave and a reflected secondary shock from the 
back wall pass through the SEe band. We apply a variety of initial perturbations to the interfaces 
between the two fluids in which the physical standard deviation, wave number range, and the spectral 
slope of the perturbations are held constant, but the number of modes initially present is varied. By 
thus decreasing the density of initial spectral modes of the interface, we find that we can achieve as 
much as 25% less total mixing at late times. This has potential direct implications for the treatment 
of initial conditions applied to material interfaces in both 3D and reduced dimensionality simulation 
models. 


I. INTRODUCTION 

Mixing of materials due to turbulent behavior has been 
a long-standing problem of interest in many areas of 
physics. In particular, the growth of interface pertur¬ 
bations due to impulsive driving of an interface between 
materials has been of great interest in diverse subjects 
including inertial confinement fusion [T], relativistic as- 
trophysical jets [2], and supernovae remnants [3]. All 
of these subjects include mixing across interfaces due to 
shock-driven turbulence. The Richtmyer-Meshkov insta¬ 
bility (RMI) [4] drives the growth of initial perturbations 
until amplitudes become sufficient to drive nonlinear in¬ 
teractions between modes. This coupling of modes drives 
a transition from laminar to turbulent behavior and thus 
allows turbulent mixing of the two fluids. RMI adds the 
challenges associated with shock physics and transitional 
turbulence to the already difficult problem of hydrody- 
namical mixing. 

Small-scale resolution requirements for simulations 
typically focus on those of continuum fluid mechanics 
described by the Navier-Stokes equations; different re¬ 
quirements are involved depending on the regime consid¬ 
ered and on the relative importance of coupled physics 
such as multi-species diffusion and combustion as de¬ 
termined by Knudsen, Reynolds, Schmidt, Damkohler, 
and other characteristic non-dimensional numbers; on 
the other hand, the longest wavelengths that can be re¬ 
solved are constrained by the size of the computational 
domain. Ideally, we would like to resolve all relevant 
space/time scales and material interfaces in our simu¬ 
lations, the so-called direct numerical simulation (DNS) 
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strategy. DNS is prohibitively expensive in the foresee¬ 
able future for most practical flows and regimes of inter¬ 
est at moderate-to-high Reynolds number (Re). On the 
other end of the simulation spectrum are the Reynolds- 
Averaged Navier-Stokes (RANS) approaches which are 
the preferred industrial standard. In coarse grained sim¬ 
ulation (CGS) strategies, large energy containing struc¬ 
tures are resolved, smaller structures are filtered out, and 
unresolved sub-grid scale (SGS) effects are modeled; this 
includes classical large eddy simulation (LES) strategies 
with explicit use of SGS models, [5] and implicit LES 
(ILES) [6], relying on SGS modeling and filtering pro¬ 
vided by physics capturing numerical algorithms. The 
CGS strategy of separating resolved and SGS physics 
effectively becomes the intermediate approach between 
DNS and RANS. 

Turbulent material mixing can be usefully character¬ 
ized by the fluid physics involved: large-scale entrain¬ 
ment, stirring due to velocity gradient fluctuations, and, 
molecular diffusion. At moderately high Re, when con¬ 
vective time-scales are much smaller than those asso¬ 
ciated with molecular diffusion, we are primarily con¬ 
cerned with the numerical simulation of the first two 
convectively-driven processes. These processes can be 
captured with sufficiently resolved ILES [6] , which serves 
as our primary simulation strategy. 

While the RMI has been studied in a wide variety of 
geometries IMO], here we will focus on a shock-tube con¬ 
figuration with a band of high density gas (sodium hex¬ 
afluoride SEe) embedded in a low density gas (air). This 
will allow us to study mixing from the initial shock, the 
rarefaction wave, and reshock after primary shock reflec¬ 
tion off the back wall of the system. 

In instability-driven turbulence, special attention has 
been paid to the initial conditions (IC) which seed the 
ballistic and nonlinear phases of instability growth and 
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the transition to turbulent behavior. Extensive numeri¬ 
cal and laboratory studies have shown that the energetic 
and mixing properties of RMI-induced flows carry signifi¬ 
cant imprint from their IC even at late times. Beyond the 
work surveyed in [4] , investigation of IC effects on RMI 
have been subject of many experimental, [nHH] com¬ 
putational, [g 0 HHH] as well as theoretical studies. 
0122]. These investigations follow on the current recog¬ 
nition that depending on the IC, different far field or late 
time self-similar regimes are possible [23]. In particular, 
computational |9| and laboratory [14] experiments show 
that the addition of high-wave number modes to the ini¬ 
tial material interface perturbations can greatly increase 
late time mixing after reshock. 

In this paper, we extend this investigation of the de¬ 
pendence of late-time behavior on IC by modeling the 
impact of a Mach number 1.26 shock on a 15 cm wide 
band of SFq. For each simulation, we apply initial per¬ 
turbations to both the front and back surfaces of the 
SFe band. These perturbations have the same physical 
standard deviation, range of wave numbers, and spec¬ 
tral slope in all simulations presented in this paper. We 
examine the effect of varying the initial spectral density 
by modifying spectral content over a given band of wave 
numbers using fractions of the total number of modes 
possible. We find that in cases with fewer initial Fourier 
modes we can show as much as 25% less mixing at late 
times compared to simulations with a full spectrum of 
initial perturbations. We attribute this variation to the 
hinderance of the development of a turbulent cascade. 


A. Guidance from Laboratory Experiments 

Laboratory studies of RMI-driven mixing provide cru¬ 
cial guidance for the simulations which we present in this 
work. A rich history of experimental work has shown that 
RMI-induced flows and mixing are strongly dependent on 
the shock velocity, the geometry of the system, and the 
initial interface perturbations M- Extensive effort has 
yielded a wide variety of measurement techniques which 
allow the extraction of rich datasets from experiments, 
including 2D concentration maps which enable measure¬ 
ments of the mixing rate [25], and simultaneous concen¬ 
tration and velocity field diagnostics for detailed analysis 
of turbulent mixing [26] . 

With advanced diagnostic tools, laboratory experi¬ 
ments have been designed to study the impact of shock 
velocity and IC on the mixing properties of RMI-driven 
turbulence. IC have been shown to play a significant role 
in the development of the flow and mixing properties of 
experiments. In shock-tube experiments, inconsistencies 
between realizations has been attributed to variations 
in IC [12] [14] . These studies provide important general 
trends against which we can compare our results. 


B. Numerical Investigations 

Simulations of RMI-induced mixing are fundamentally 
limited by the nature of shock discontinuities. Without 
the possibility of direct numerical simulations (DNS), 3D 
models are limited to a large-eddy simulation (LES) ap¬ 
proach in which large scales are computationally resolved 
and the effects of unresolved scales is modeled either ex¬ 
plicitly or by the choice of highly-stable numerical oper¬ 
ators. In somewhat similar fashion, laboratory experi¬ 
ments are limited by the finite scales of the instrumental 
resolution of measuring/visualizing devices characteriz¬ 
ing the resolution of the observational techniques EZj. 

In the LES framework, unresolved scales are generally 
assigned a diffusive character and represented by a suit¬ 
able diffusive operator. One option is to add explicit 
diffusive terms which represent unresolved diffusive pro¬ 
cesses. A wide range of such explicit operators have been 
studied in a variety of engineering, geophysical, and as- 
trophysical settings [28H32] . 

Classical LES is particularly inadequate for flows 
driven by RMI because of the dissipative numerics needed 
for shock capturing [5]. Hybrid methods switching be¬ 
tween shock-capturing schemes and conventional LES de¬ 
pending on the local flow conditions have been proposed 
na. Such high-order shock-capturing (e.g., fifth-order 
WENO) methods are typically chosen to ensure a smooth 
transition and matching of the inherently different sim¬ 
ulation models. However, all shock capturing methods 
degrade to first-order in the vicinity of shocks because of 
the monotonicity requirements and in particular, at the 
very important initial stage at which a shock first inter¬ 
acts with a compositional interface and generate the fluc¬ 
tuating velocity field. Thus, severe resolution demands 
to address the competition between the implicit subgrid 
models provided by numerical operators and explicit sub¬ 
grid models can be expected in the hybrid context. Al¬ 
ternatively, many codes now employ so-called physics¬ 
capturing numerical schemes which by design rely on 
minimal numerical diffusion for computational stability. 
This ILFS framework permits simulations to run with 
minimal levels of diffusion which is crucial for correctly 
modeling shock physics [6] . By combining shock and tur¬ 
bulence emulation capabilities based on a single physics¬ 
capturing numerical model, ILFS provides a natural ef¬ 
fective simulation strategy for RMI [33]. Moreover, ILFS 
can accurately capture the stirring-driven scalar mixing 
which is the dominant aspect for high Re flows [34| . 

Extensive numerical studies have investigated many 
properties of the RMI in a variety of geometries, in¬ 
cluding the single planar interface 00[TSH1I1I3S|, the 
two-interface inverse chevron [Hiini, the shocked heavy 
cylinder [36], and the shocked gas curtain [37]. Spe¬ 
cial attention has been paid to the effects of IC specifics 
[9l|35]. ILFS codes have been also used in a wide vari¬ 
ety of hydrodynamic and magnetohydrodynamic applica¬ 
tions including airfoil design, planetary magnetospheres, 
inertial confinement fusion, and stellar dynamo models. 
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FIG. 1. (Color online) Schematic diagram of the IC for all simulations (see Table |^. The shock wave initially propagates to 
the right. Solid lines represent solid boundaries at both ends of the shock-tube, while dotted lines show periodic boundaries in 
the transverse directions. Initial perturbations are applied to the front and back interfaces of the SFe band, labeled Li and L 2 
respectively. 


e.g.,[TJ|38H40]. 


II. COMPUTATIONAL TOOLS AND 
SIMULATION PARAMETERS 

Here we use the ILES capability of the code RAGE 
jlT] to study the RMI and associated induced turbulent 
mixing. RAGE has been used extensively for similar sim¬ 
ulations, reporting good agreement with available labora¬ 
tory observations P m [lOl |37l |42] . RAGE solves the Eu¬ 
ler equations of multi-fluid compressible hydrodynamics 
using a secondorder Godunov finite-volume scheme, and 
offering a variety of numerical options. In the present 
simulations, adaptive mesh refinement and material in¬ 
terface RAGE options were not used; a Van Leer limiter 
and uniform Cartesian gridding were chosen. 

The evolution equations encode conservation of mass, 
material concentrations, momentum, and energy given 
by 


| + v.(,„-) = o 

(1) 

^ + V ■ ipY^u) = 0 

(2) 

^ + (w-V) ipu)+VP = 0 

(3) 

^+V-[{pE + P)u]=0, 

(4) 

where p is the mass density, u is the velocity, W is the 
mass fraction of species n, P is the fluid pressure, and E 
is the total specific energy. These equations are supple¬ 
mented with SESAME tabular equations of state [43] . 

We define the transverse average of a quantity 


f{x,y,z) as 


f{x) 


11 fdydz 
LyLz 


(5) 


and the mass-weighted transverse mean as 


hx) = 


11 Pfdydz 
pLyLz 


(6) 


We have neglected body forces, molecular shear, and heat 
conduction. The ratio of specific heats for mixtures is 
computed using a volume-fraction weighted average of 
the ratio of specific heats for each gas. As used in the 
present work, RAGE models miscible material interfaces 
(Schmidt number S'c ^ 1) and high Re convection-driven 
flow with effective numerical viscosity determined by the 
small-scale cutoff. Timesteps in RAGE runs are set as 
minimum of a Courant condition on the hydrodynamics 
and stability criteria determined by other physics in the 
code. 

Here we report on a series of nine simulations which all 
follow the basic physical setup shown in Eigure The 
simulation domain is 54 cm long by 7.62 cm square in 
the transverse direction. At the start of the simulation 
shocked air fills the leftmost 15 cm of the domain, with a 

5 cm gap between the shock front and the first interface 
between the unshocked air and the SEe, which we label 
Li. The initially stationary SFq band has a width of 
15 cm. The right side of the band is labeled the second 
interface between unshocked air and SEe or L 2 . The 
reminder of the domain is filled with unshocked air. The 
Atwood number for both interfaces is 0.69. 

A shocked air region is created upstream in terms of 
a higher-density higher-pressure region chosen to satisfy 
the Rankine-Hugoniot relations for a Mach 1.26 shock. 
The primary shock propagates in the x-direction through 
the contact discontinuities and reflects at the end of the 
simulation box on the right. Boundary conditions in¬ 
volve reflecting walls on the right in the x-direction and 
periodicity in the transverse directions. 

Eor this system, a simple estimate of the Reynolds 
number for the physical system can be found using the 
shock velocity, the transverse width of the domain, and 
the viscosity of air. This gives a Re of approximately 

6 X 10^. The enormous range of scales which would be 
needed make a true DNS calculation impractical for the 
foreseeable future. Thus we choose a CGS strategy - 
RAGE based ILES. 
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(c)Initial Surface Perturbations for Case S-8 
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(e)Initial Surface Perturbations for Case S-32 
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FIG. 2. (Color online) Initial surface perturbations with (a) 100 modes, (b) 800 modes, and (c) 3263 modes, along with the 
spectral content of those surfaces shown in (b), (d), and (f) for each case respectively. All surfaces have a standard deviation 
of 0.1 cm and random phasing of the modes. For (b), (d), and (f), red (medium gray) indicates modes with wave numbers 
that do not satisfy Equation]^ green (light gray) indicates eligible modes that are not used in this realization of the interface 
surface, and blue shows (dark gray) modes which have been used to construct the associated surface. Note that in (f) all modes 
compliant with Equation are used. 
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TABLE L Simulation parameters for primary models, includ¬ 
ing the numbers of grid points in the x (streamwise), y, and 
2 : directions Nx x Ny x Nz, the grid spacing Ax in /xm, and 
the number of modes used to construct the initial interface 
perturbations on either side of the SFe band, and the number 
of randomly generated realizations simulated Nr. 


Case 

Nx X Ny X Nz 

Ax (/xm) 

Nic 

Nr 

V-1 

4552 X 600 X 600 

127 

100 

1 

H-32 

2990 X 422 x 422 

181 

3263 

1 

H-1 

2990 X 422 x 422 

181 

100 

1 

S-32 

1800 X 254 X 254 

300 

3263 

3 

S-16 

1800 X 254 X 254 

300 

1600 

1 

S-8 

1800 X 254 X 254 

300 

800 

1 

S-4 

1800 X 254 X 254 

300 

400 

3 

S-2 

1800 X 254 X 254 

300 

200 

1 

S-1 

1800 X 254 X 254 

300 

100 
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A. Initial Interface Perturbations 

IC are applied to both interfaces as perturbations in 
the X position of the interface at each point. Thus the 
perturbed interface location is given by where 

is the perturbation of the interface. We define 7 /^ 
as 

Nic 

V’ {y, ^) = 

+ A sm{kiy + kz)] 

(7) 

where Njc is the number of modes to be applied as initial 
perturbations, Ai is the amplitude of mode with wave 
number ki in the ^-direction and li in the ^-direction. 
Coefficients Ci and Di are randomly selected from a 
normal distribution with 0 mean and a variance of 1. 
This is equivalent to stating that the modes have purely 
random phases. Amplitudes Ai are chosen such that 
Ai oc {k^+nr so that the amplitude of a given mode 
scales as the magnitude of its wavevector to the —2 
power. The values for all Ai are uniformly scaled so that 
the standard deviation of the interface perturbations is 
0.1 cm for all simulations. Different mode sets with iden¬ 
tical standard deviation, spectral slope, and numbers of 
initial modes are used for the the front and back inter¬ 
faces of the band to avoid any type of synchronization 
between the two interfaces. 

Our broadband initial perturbations all utilize the 
same range of wave numbers. We select modes such that 

3.298 cm-i < + < 52.77 cm-^ (8) 

or equivalently 

4 < ml Nnl < 64 (9) 



FIG. 3. (Color online) A time-space diagram showing the 
position of the front (bottom solid line) and back (second solid 
line from bottom) of both mixing layers Li and L 2 (front is 
second solid line from top, back is top solid line) as well as 
the primary shock and the rarefaction wave for case S-32. 


where and rii are the mode numbers in the y- and z- 
directions, respectively. This means that there are 3263 
possible distinct wave number pairs that can generate 
orthogonal modes, which we can then apply as interface 
perturbations. For simulations with fewer than the max¬ 
imum number of modes, we randomly select a number of 
modes Njc from the 3263 total possible. The selection 
process is done by choosing Nic random integers from 
a uniform distribution between 0 and 64 for and rii 
independently such that the constraint in Equation is 
satisfied. We also do not allow duplicate modes, so acci¬ 
dentally occurring repeats of {mi, rii) pairs are removed 
and reselected. In this way the modes correctly sam¬ 
ple the mode number space with the same fraction of 
possible modes selected on average from any two bands 
in mode magnitude. For example, the interface per- 
turbatio ns used with Njc = 100 had 10 modes with 
10 < ml + nl < 20 out of a possible 249 modes, or 
3.7% of all possible modes, and 15 out of 401 modes with 
mode magnitudes between 20 and 30, or 4.0%. As ex¬ 
pected, cases with Njc = 200 utilize 6.8% and 7.0% 
of the possible modes in these two bands, respectively. 
The percentages of utilized modes are 12.9% and 13.2% 
for cases with Njc = 400, 25.3% and 24.2% for cases 
with Nic = 800, and 49.8% and 49.4% for cases with 
Nic = 1600. Thus the fraction of possible modes utilized 
in any band in mode magnitude will be roughly equal to 
the overall fraction of possible modes used. We define 
this fraction of modes used over any range of values in 
mode magnitude as the spectral density of the surface. 

Figure shows the initial interface perturbations ap¬ 
plied to the front interface Li for simulations with 100, 
800, and 3263 modes, as well as a schematic representa¬ 
tion of the modes present in each interface. All interfaces 
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FIG. 4. (Color online) Snapshot 3D volume renderings of SFe mass fraction from x = 18 cm to x = 54 cm for case S-32 at (a) 
the start of the simulation, (b) t = 1.0 ms when the shock front is in the SFe band, (c) t = 2.0 ms when the shock has exited 
the band, and (d) t = 3.0 ms after reshock. Ysfg is shown with low concentrations in yellow (light gray), high concentrations 
in purple (dark gray), and pure air transparent. Visualizations made with VAPOR [44] . 


have the same RMS amplitude, but visual differences can 
be noticed due to the differing spectral content. For ex¬ 
ample, comparing Figures [^a) and (e) reveals that the 
case with 32 times the number of modes shows consid¬ 
erably more spatial inhomogeneity. Quantitatively, the 
surface with Njc = 100 has a kurtosis of 2.24 while 
the surface with Njc = 3263 has a kurtosis of 3.02, or 
slightly more than Guassian. As the number of modes in¬ 
creases, the surface perturbations become more Gaussian 
and thus approach a kurtosis of three. 

Table |T] list the simulations discussed in this paper and 
their attributes. For all but cases V-1, H-32 and H-1, 
we use a grid spacing of 0.03 cm, giving a computational 
domain that is 1800 cells by 254 cells square for a total 
of 116 million computational elements in each simula¬ 
tion. Cases H-32 and H-1 use a grid spacing of 0.018 cm 
with the same physical domain for a total of 929 million 
computational elements. Case V-1 uses a grid spacing of 
0.0127 cm for a total of 1.64 billion computational cells. 

As both the modes selected and the phases of each indi¬ 
vidual mode are randomly assigned there is the potential 
for different realizations of our IC to produce significantly 
different results. We might expect this realization noise 
to be particularly pronounced for low numbers of initial 


modes where there may be only a few modes selected at 
the longest wavelengths. To address this issue, we have 
run multiple realizations of cases S-32, S-4, and S-1 in 
order to give some idea of the magnitude of the realiza¬ 
tion noise at these three values of Njc discussed further 
below. 


B. Simulation Results 

All simulations were begun from IC like those illus¬ 
trated in Figure and evolved for at least 3.0 ms. Fig¬ 
ure shows the evolution of case S-32 as a sample. The 
primary shock hits the first interface about 0.1 ms after 
the start of the simulation and begins compressing the 
band. The mixing layer at this interface begins to grow 
at that point. The primary shock hits the back inter¬ 
face at about t = 0.92 ms, launching a rarefaction wave 
and beginning the growth of the second mixing layer. 
The primary shock reflects off the back wall at about 
t = 1.38 ms. The rarefaction wave hits the first interface 
at about t = 1.50 ms and continues to propagate through 
the shocked air. The reflected primary shock hits the sec¬ 
ond interface at about 1.55 ms and proceeds to reshock 
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the band until it hits the first interface again at about 
t = 2.15 ms. 

Following previous studies of shock-driven mixing [9l 
[M [371 [45], we can describe the mixing layers at each 
interface in terms of the SFq mass fraction, Isfg, using 
the mixing parameter 0 defined as 

= 41sf 6 (1 - hsFe), (10) 

where (p ranges from zero for pure air or SFe to one for 
a 50/50 mix by mass fraction of both species. In dis¬ 
cussing material mixing in shock-driven turbulence, there 
are a wide variety of metrics used to measure mix growth. 
Inspired by reduced dimensionality RANS models, one 
commonly used metric is the width of the mixing layer 
or mix width. Here we use the product transverse aver¬ 
ages of the mass fractions of air and SFq to define the 
mean mixing p{x) as 

0(rc) = (1 - . (11) 

We define the mix width W as the length of the interval 
in the stream-wise direction where p{x) > 0.75. While 
this definition is somewhat arbitrary, we have found that 
our results do not vary significantly when values of 0.7 
or 0.8 are used. Our choice of 0.75 has been exten¬ 
sively used before [9] and is designed to focus the anal¬ 
ysis on the regions of strong mixing while also providing 
wide enough layers to limit realization noise for volume- 
averaged quantities. 

Using this metric, we can define the volumes of pri¬ 
mary mixing for both interfaces. We label the mixing 
layer associated with the first or left interface Li and the 
mixing layer associated with the second or right inter¬ 
face 1/2- In Figure the two mixing layers are shown 
to evolve in width over time. We will make particularly 
detailed study of these two layers in subsequent analysis. 

Figure shows 3D volume renderings the SFe concen¬ 
tration at t = 0.0, I.O, 2.0 and 3.0 ms for case S-32, high¬ 
lighting the movement of the band as well as the growth 
of both interfaces. Of particular note is the transition 
from ballistic growth seen at t = 2.0 ms to the strong 
mixing seen by t = 3.0 ms. While the quantitative de¬ 
tails of these simulations differ, all follow the same basic 
evolution outlined here. 


III. RESOLUTION STUDY 

In this section, we compare cases H-32 and S-I, and 
V-I, H-I, and S-I, respectively, to address convergence 
issues. Firstly, we note that previous simulations using 
RAGE for a single planar-interface geometry have shown 
that large scale convergence to within chaotic variabil¬ 
ity can be achieved provided the initial interface per¬ 
turbations are well resolved [9]. Robustness of ILFS 
macroscopic results is expected with sufficient resolution 
once there is sufficient separation in scales between the 
energy-carrying scales and the diffusive scales in the high 


Reynolds number limit [34| . To investigate the degree to 
which our high-Re assumption is true, we have conducted 
several simulations at increased resolution. 

Secondly, it is possible that the effects of varying the 
initial spectral density may change when the dissipation 
of the system is changed. In ILFS codes such as RAGE 
the effective (numerical) diffusion depends on the grid 
resolution. [34l [46] Thus our resolution study is actually 
a study of both increased resolution and decreased effec¬ 
tive diffusion. In a well-developed turbulent system the 
dissipation seen on moderate scales should become inde¬ 
pendent of the scale at which dissipation occurs, as the 
cascade reaches a steady-state of transfer. However, for 
transitional flows this is less certain. In running cases 
V-I, H-I, and S-I, and again in cases H-32 and S-32, 
identical initial conditions were used for both interfaces. 
Unlike the independent realizations of our initial condi¬ 
tions that were run for cases S-I, S-4, and S-32 which 
used unique random spectra for each interface in each 
realization, here we have used identical modes with iden¬ 
tical phases in order to isolate the effects of resolution 
in our simulations. Ideally it would be desirable to also 
run multiple realizations (with randomly varying mode 
weights) at each resolution, however the computational 
demands of H- and V-series cases would make such vari¬ 
ability studies impractical within our present scope; how¬ 
ever expected sensitivity of results to such variations is 
demonstrated further below. 


A. Reynolds Number and Spectral Diagnostics 

When studying the resolution dependence of our re¬ 
sults, we chose to first examine the quantities we would 
expect to vary in our models, namely Re and the kinetic 
energy spectrum. Only estimates for Re for these sim¬ 
ulations are possible as our ILFS framework does not 
include physical dissipation. Furthermore any discussion 
of Re must account for the transitional nature of the 
flows as the RMI develops from ballistic growth to mode 
interactions to turbulence and finally to decaying turbu¬ 
lence. We can estimate an ILFS effective Re by finding 
a measure of the effective viscosity z^eft associated with 
resolution estimated as ratio of computed dissipation and 
squared strain-rate magnitude using the raw simulated 
flow velocity [34] - distinct from the grid-scale viscosity 
associated with the numerical scheme specifics. Zhou et 
al. [46] have explored several methods for estimating the 
Here we will utilize two methods, which both posit 
that the effective viscosity can be represented by 

t'eff = ^, (12) 

where e is the energy dissipation rate and Q, = l/2(Vxi7)^ 
is the enstrophy. The first method, which we will use to 
calculate what we term Rei uses the resolved energy flux 
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FIG. 5. (Color online) (a) Ensemble-averged Re for the S- and H-series cases, along with the Re number for case V-1. The Re 
for all three resolutions show significant variability in time, however there is a clear increase at late times for higher resolution 
simulations, (b) Streamwise-averaged kinetic energy power spectra for cases S-1 (red/medium gray), H-1 (green/light gray), 
and V-1 (blue/ dark gray). All spectra are computed for wave numbers less than 16.7 cm“^. All cases show similar spectral 
features at low wave numbers, but increasing resolution leads to extend inertial ranges. 


spectrum 


e 


= 




rk* dE{k,t) 

Jo dt 




dk 


(13) 


where /c* is the wave number chosen at which the energy 
transfer rate has achieved a constant value, and E{k^ t) is 
the Fourier transform of the specific kinetic energy [46]. 
Using this prescription, we can then compute Rci as 


Rei 


um 


rk* dE{k,t) 

Jo dt 


dk 


(14) 


where U is the volume-averaged magnitude of the fluctu¬ 
ating velocity over the mixing layer, L is the width of the 
mixing layer, and U is the volume-averaged enstrophy. 

The second method explored by Zhou et al. m uses 
the dimensionless parameter which is given by 



(15) 


thus providing an alternative means of computing the 
dissipation rate e in terms of the characteristic velocity 
U and length L scales if D is known [46]. DNS models 
of forced and decaying turbulence have shown that D 
asymptotically tends to a constant value of order unity 
with increasing Re [4^ . Assuming our Re is high enough 
that our simulations fall in the asymptotic regime, we 
can compute another estimate for the Re as 


Re2 


IP' 


(16) 


Here, we have again used the mix width for our length 
scale L, the volume-averaged magnitude of the velocity 
fluctuations as our velocity scale U, and the volume- 
averaged enstrophy U. 

Interestingly, we find that both estimates provide very 
similar values for each simulation as a function of time. 
We compute what we term a mean Reynolds number Re 
as 


Re — — (Rci + Re 2 ). (17) 

We further compute ensemble means for all cases at a 
given resolution at each time in order to provide an ad¬ 
ditional reduction in noise. Figure |^a) shows the mean 
Re for each resolution as a function of time. While there 
is considerable noise in these calculations, there is also a 
clear trend for higher resolutions to produce larger Re at 
late times, approaching the (integral-scale based Re) mix¬ 
ing transition threshold values between 10^ and 2 x 10^ 

m- 

An additional diagnostic of our resolution study is ex¬ 
amining the length of the simulated inertial range in the 
kinetic energy power spectra of our mixing layers. Higher 
Re flows should exhibit a greater separation in scales be¬ 
tween the energy-containing large scales and the dissi¬ 
pation range, and we should find a larger self-similar 
inertial subrange with expected Kolmogorov scaling of 
kinetic energy power law ~ k~^P. Figure |^b) shows 
the fluctuating kinetic energy power spectra for the back 
mixing layer at t = 2.8 ms. Case V-1 clearly shows the 
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FIG. 6. (Color online) (a) Turbulent kinetic energy (TKE) as a function of time for simulations S-32, S-1, H-32, H-1, and V-1. 
(b) Total mixednes (TMX) for the same simulations. Note that lines for all five cases in both plots are present at all times, 
though they are often over-plotted by another case as values are very similar. For both quantities the evolution does not show 
significant changes as the numerical resolution is increased and the variation caused by changing the initial spectral density is 
considerably larger than the variation due to the change in resolution. In particular, cases H-1 and V-1 show less than 2.8% 
variation in TKE and less than 1.1% variation in TMX at late times. 


longest inertial range, while case H-1 and case S-1 suggest 
progressively smaller inertial ranges. 


B. Large-Scale Behaviors 


While Re and small-scale behaviors of our models 
clearly do vary with resolution, we will show that the 
S-series cases are sufficiently resolved in to achieve some 
level of convergence on the large-scale features of interest 
to us as we increase resolution in cases H-1, H-32, and 
V-1. Specifically, we will examine the volume-integrated 
turbulent kinetic energy and mixedness as a function of 
resolution. 

Figurej^a) shows the evolution of the turbulent kinetic 
energy TKE over the course of the simulation for five 
cases. We define TKE as 


TKE = 



- 


+ («!/ - '^vf + 



dV, 


(18) 

where V represent the entire simulated volume. TKE 
thus removes the kinetic energy associated with any flows 
along the axis of the shock tube, most importantly the 
shocks themselves. We have plotted TKE for cases V-1, 
H-32, H-1, S-32, and S-1. Cases H-32 and S-32 show very 
similar overall evolution with some deviation, as would 
be expected for chaotic systems like ours. Cases V-1, H-1, 
and S-1 also show very similar overall behavior. Careful 
examination does show that for both values of Nic the 


S-series resolution simulations both consistently show be¬ 
tween 2% and 5% more mixing after reshock than their 
H-series counterparts, indicating that there may be some 
aspects of this problem that are not fully numerically con¬ 
verged, however when compared to the variation caused 
by changes to the initial spectral density the variations 
due to increased numerical resolution are between 3 and 
6 times smaller. Case V-1 shows less than 3% variation 
from case H-1 at all times and no clear trend relative to 
case H-1. We attribute this variation, which is largest 
after 1.6 ms, to chaotic variability. 

Figure ib) shows the evolution of the volume- 
integrated mixedness of the same five simulations. In 
our RAGE simulations, the composition inside a compu¬ 
tational cell is assumed to be uniform throughout that 
cell, thus materials are numerically treated as if they are 
a uniform molecular mix. We term this mixedness. One 
measure of mixedness is the volume-integrated, density- 
weighted mixedness TMX. Like TKE, TMX is a volume- 
integrated measure of the mixing between the air and 
SF6 in our domain. TMX is defined as 

TMX= [ p^YsmY^^dV, (19) 

Jv 

where Fair can also be expressed as 1—1 sf 6- Thus TMX is 
the product of the partial mass fractions for the two ma¬ 
terials in a computational cell at any instant, integrated 
over the computational domain. We can consider the 
volume integrated TMX over the whole domain, thereby 
including both interfaces in the same metric, or we can 
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(b) Relative TKE 


FIG. 7. (Color online) (a) Turbulent kinetic energy as a function of time for simulations S-32 through S-1. (b) Relative TKE 
for the same simulations where the ensemble mean at each time from (a) has been subtracted for each simulation. While strong 
variations are seen, there is no obvious trend with Nic- 


divide the domain into two halves at the center of the 
SFe band, thereby allowing us to examine the mixed¬ 
ness for each interface separately. We also examine the 
normalized total mixedness, where the normalization is 
provided by the TMX computed for our IC. 

In an ILES code such as RAGE, increases in mixedness 
are only facilitated by numerical diffusion of species con¬ 
centration gradients at or near the grid scale. The numer¬ 
ical operators in RAGE are essentially inviscid in regions 
without small-scale gradients. In our IC, only cells on 
the interfaces at the front and back of the SFe band have 
any mixedness. All subsequent growth in TMX is thus a 
useful diagnostic of the effectiveness of the turbulent cas¬ 
cade in moving energy (and hence material interfaces) to 
sufficiently small scales where numerical dissipation can 
act. Generally this numerical diffusion is minimized in 
an ILES framework, but is still likely to be much greater 
than realistic molecular diffusion for this application. 

As with TKE, Figure [^b) demonstrates that TMX 
for all five cases follow the same general trend with slow 
growth prior to t = 1.6 ms followed by much more rapid 
mixing. Cases H-32 and S-32 show very similar trends 
with deviations of less than 4% over their entire evolu¬ 
tion. Cases H-1 and S-1 also show highly similar evo¬ 
lution with maximum deviation of only 6%, while cases 
H-32 compared with H-1, and S-32 compared with S-1 
show much greater disagreement of as much as 18% and 
16%, respectively. In the total mixedness, the variation 
between cases H-1 and V-1 is very small. At t = 2.0 
ms case V-1 shows 0.73% less mixedness than case H-1, 
while at t = 3.0 ms case V-l’s TMX is 1.12% smaller 
than case H-l’s. Unlike TKE, TMX does show a very 
small but systematic decrease between cases H-1 and V-1. 


We speculate that continued increases in resolution may 
yield asymptotically smaller reductions in TMX. Thus for 
the purposes of examining the impact of initial spectral 
density on mixing, we feel confident that our numerical 
resolution is sufficient. 

For both global measures, increased resolution has 
smaller changes to the volume-integrated quantities mea¬ 
suring turbulent kinetic energy and mixedness. Thus we 
believe that the primary effect in our models is the change 
in initial spectral density rather than any issues arising 
from insufficient numerical resolution. 


IV. EFFECTS OF SPECTRAL DENSITY OF 
INITIAL PERTURBATIONS 

The primary purpose of this paper is to investigate 
the effects of the spectral content of initial interface per¬ 
turbations on the growth of the RMI, the coupling of 
modes through nonlinear behaviors, and the development 
of turbulence. To that end, we have conducted a series 
of otherwise identical simulations with variable numbers 
of Fourier modes in their initial interface perturbations. 
These simulations are labeled S-1, S-2, S-4, S-8, S-16, and 
S-32 (see Table . All six cases follow the same general 
time evolution as shown for case S-32 in Figure how¬ 
ever significant variations are seen in the properties of 
the resulting mixing layers around both interfaces both 
before and after reshock. 

Figure [^a) shows the evolution of the total turbulent 
kinetic energy (TKE) in the computational domain. Ex¬ 
amining the time-evolution of TKE, we can see sharp 
increases at 0.1 ms when the shock first impacts Li, at 
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FIG. 8. (Color online) (a) Time-evolution of the width of the mixing layer Li for cases S-1 through S-32. The mix widths 
grow slowly until the rarefaction wave hits it at about t = 1.6 ms, after which all six cases show a rapid expansion of the 
layer. This expansion is reversed upon reshock, which recompresses the the layer, (b) Time-evolution of the mix width for the 
second interface L 2 for all six cases. The mix width grows slowly until reshock at about t — 1.6 ms where a sharp compression 
occurs and is followed by continued expansion, (c) Relative mix width for the front mixing layer Li where the ensemble mean 
of the six cases at each time has been subtracted. No clear trends with Nic can be seen either in the orderly evolution prior 
to reshock or the more chaotic behavior after reshock, (d) Relative mix width for the back interface L 2 where the ensemble 
mean of the six cases has again been subtracted at each time. While the mix widths are less chaotic than in Li, there is no 
trend with Nic- 


0.9 ms when the main shock impacts I/ 2 , around 1.6 ms 
when the rarefaction wave hits Li, and about 1.6 and 2.2 
ms when the reflected shock excites additional motion at 
each interface. All six cases follow the same general tim¬ 
ing, however the relative levels of TKE achieved differ 
widely. 

Figure mb) shows evolution of TKE in all six cases 
with the ensemble mean at each time subtracted. This 


highlights the variation is the six cases. Variations of as 
much as 50% can be seen between cases, most notably 
between cases S-4 in brown and S-2 in red. This is strik¬ 
ing as all cases have IC with equal spectral slopes and 
surface standard deviations. The only difference here is 
the number of modes applied as interface perturbations. 
There does not, however, appear to be a systematic trend 
with Njc at either early or late times. 
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(c)TMX by Interface 

FIG. 9. (Color online) (a) Volume-integrated mixedness TMX for cases S-32 through S-1 as a function of simulation time. 
Values have been normalized by the values of TMX for each case just prior to the first impact of the shock on the front interface, 
which we term TMXq. Two clear phases of mixing can be seen with mixing before 1.5 ms showing linear behavior followed by 
much more rapid, turbulent mixing, (b) Relative TMX for the same simulations with the ensemble mean subtracted at each 
time, (c) TMX for each interface for cases S-1, S-4, and S-32. While the initial linear phase shows little dependence on V/c, 
the later turbulent phase reveals a clear trend where cases with greater Nic exhibit more mixing. By t = 2.5 ms there is a 
monotonic relation between Nic and TMX for both interfaces individual and over the entire domain. 


As discussed above, the mix width W is defined as the 
length of the interval in the stream-wise direction where 
^{x) > 0.75. Figure shows the widths of the mixing 
layers at both the front and back of the SFe band as a 
function of simulation time. Figure [^a) shows the extent 
of Li, while panel (b) shows the width of I/ 2 . For both 
mix widths, we see slow growth after the primary shock. 


For Li we see a sharp increase in the mix width upon rar¬ 
efaction and then a sharp decrease as reshock compresses 
the mixing layer. For L 2 we see a brief compression of the 
mixing layer on reshock followed by continued expansion 
of the mixing layer. All six cases show similar behaviors 
on both interfaces. 

Figure l^c) shows the relative mix width for Li, while 
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FIG. 10. (Color online) Volume-integrated mixedness at t = 
3.0 ms for six values of Nic- Unique randomly generated 
instances of the initial perturbation spectra were used for each 
simulation, with five realizations of case S-1 and three each 
for cases S-4 and S-32. For case S-32, all three symbols are 
plotted here though they are visually difficult to distinguish. 
While there is some evidence for significant variability due 
to realization noise, particularly for small values of V/c, the 
larger trend of less mixedness with decreasing Nic is clear. 


panel (d) shows the same for L 2 . To construct the rela¬ 
tive mix width we subtract the ensemble mean for all six 
cases at each time from each curve, thus removing the 
general trend to focus on the variations between cases. 
Figure [^c) shows that while there are significant vari¬ 
ations in the mix widths of as much as 20% prior to 
rarefaction, there are no clear trends as the initial spec¬ 
tral density is changed. Particularly after rarefaction and 
reshock, the mix width of Li shows dramatic and chaotic 
variability between the six cases. The relative mix widths 
for 1/2 are somewhat less chaotic, but still no trend is seen 
with initial spectral density. 

Our turbulent mixing measures are usefully character¬ 
ized by the two different mixing processes captured with 
ILES. Mix width is largely due to large-scale entrain¬ 
ment promoted by bubbles and spikes, while mixedness 
is driven by stirring associated with smaller scale veloc¬ 
ity gradient fluctuations. In Figure [^a) we can see that, 
much as with TKE, the time evolution of TMX follows 
the same general trend for all six cases. There is a slowly- 
growing phase up until about 1.6 ms when the first rar¬ 
efaction wave hits the front mixing layer leading to the 
large jump in TKE seen in Eigurej^a) at that time. After 
this point TMX grows much more rapidly for all cases. 

The variations in TMX between the six cases are high¬ 
lighted in Eigure |^b) where we have removed the en¬ 
semble mean at each time, thus eliminating the general 


trend in TMX with time. There is a slow divergence be¬ 
tween the six cases until t = 1.6 ms when reshock of L 2 
and rarefaction of Li occur. After t = 1.6 ms the diver¬ 
gence between the cases accelerates. Variations between 
cases of as much as ±30% from the ensemble mean are 
seen even at early times. TMX shows a strong depen¬ 
dence on Nic times. By t = 3.0 ms, the relation 

between Nic and TMX is monotonic. The more dense 
the spectrum of initial perturbations applied to the inter¬ 
faces, the more mixing at late times. Eigure |^c) shows 
that this trend holds for both interfaces individually. The 
mixedness of both the front and back of the band see lit¬ 
tle or no dependence on Njc until t = 16.0 ms when 
the rarefaction wave hits the front interface and the back 
interface is reshocked. Mixing at the front interface is 
again accelerated when it too is reshocked. 

Case S-8 is particularly interesting as it shows the least 
TMX of all six simulations at t = 1.6 ms but rapidly 
increases to pass cases S-1, S-2, and S-4 hy t = 2.4 ms. 
This indicates that while mixing in the early growth of 
the interface perturbations due to the RMI might not 
be as strongly dependent on the spectral density of the 
initial perturbations as the late stages after rarefaction 
and reshock when the flows become much more turbulent. 

The dependence of TMX at late times on Njc shown in 
Eigure naturally brings up the question of realization 
noise due to the random selection of modes and their 
phases for our IC. To help address this issue we have run 
additional random realizations of cases S-32, S-4, and S- 
1. Eigure pT| shows TMX at t = 3.0 ms for the S-series of 
simulations, including five unique, randomly generated 
realizations of case S-1 and three realizations each for 
cases S-4 and S-32. As a liberal estimate of the realization 
noise for each value of Njc, we take the largest difference 
between TMX fora given value of Njc and divide by the 
average over all realizations. This yields realization noise 
levels of 0.2% for case S-32, 0.8% for case S-4, and 6.7% 
for case S-1. It is not surprising that realization noise 
scales inversely with Njc as the cases with fewer modes 
are more susceptible to changes in any single mode. This 
suggests that our reported relationship between initial 
spectral density and late-time mixedness is not simply 
an artifact of particular realizations of our random initial 
perturbations and phases. 

Another diagnostic of the mixing present in these simu¬ 
lation is the normalized probability distribution function 
(PDE) of the mass fraction of SFq in the mixing layers 
at the front and back of the band. We define these mix¬ 
ing layers using the mixing parameter (j) defined in Equa- 
tionflOl We can further examine the evolution of the PDE 
for each simulation over time. Eigure pT] shows the PDEs 
for lsF6 over the evolution of cases S-1 and S-32. Eor sim¬ 
plicity we have chosen to present only the two limiting 
cases, but the intermediate simulations show behaviors 
which both qualitatively and quantitatively fall between 
there two cases. Eigure pT]^a) shows the time evolution 
of the PDE for 1 sf 6 for case S-32 with 30 bins at 0.1 
ms intervals. Lighter colors indicate larger fractions of 
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FIG. 11. (Color online) Time-evolution of the normalized probability distribution for Ysf6 in mixing layer Li for (a) Case S-32 
and (b) Case S-1 with light tones indicating larger fractions of the mixing layer at a given concentration and dark tones showing 
smaller fractions of the mixing layer volume, (c) The ratio of (b) to (a), with blue to green (darker gray) tones indicating larger 
values for case S-32 and red to yellow (lighter gray) tones indicating larger values for case S-1. 


the mixing layer Li. Large fractions of the mixing layer 
are dominated by either pure SFe or pure air until 2.5 
ms when the most of the mixing layer rapidly becomes 
mixed with 1 sf 6 between 0.4 and 0.8. Figure [T^b) re¬ 
veals qualitatively similar behavior for case S-1. 

Figure [T^c) shows the ratio of the PDFs over time in 
Case S-1 to those in Case S-32, which blue tones indicat¬ 
ing larger values for a given bin at a given time in Case 
S-1 and red values showing larger values for Case S-32. 
Prior to 2.5 ms, the primary difference between these two 
simulations is that Case S-1 shows somewhat more vol¬ 
ume with lsF6 ^0.6 while Case S-32 shows more volume 
for lsF6 ^0.6. After 2.5 ms, however, sharp differences 
become more apparent. Case S-1 shows a strong excess 
of volume with pure or nearly-pure materials compared 
to Case S-32, with more than 50% larger fractions of 
the mixing layer in the lowest and highest bins. Case 
S-32 shows as much as 50% larger fractions of the mixing 


layer with 0.1 ^ 1 ^sf6 ^ 0-4 and roughly equal values for 
4 sf 6 ^0.5. From this we can see that Case S-1 includes 
more “blob”-like behavior of the two fluids as resolved 
nearly-pure chunks of both gases pass each other, while 
Case S-32 shows much more mixing of the two gases at 
the grid-scale. 

Figure pr] shows that increased initial spectral density 
leads to enhanced mixing at the grid-scale. This may be 
due to the inhibition of a turbulent cascade of energy to 
small scales in cases where the initial spectral density is 
low. In effect, the enhanced mixing may be the result 
of the IC in case S-32 more closely approximating the 
cascade established naturally by turbulence. To test this 
concept, we need to explore the spectral content of the 
irregularities in the mixing layers. 

We construct the streamwise-averaged transverse ki¬ 
netic energy power spectra by first isolating the volumes 
for each mixing layer, then taking the 2D Fourier trans- 
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FIG. 12. (Color online) (a) Streamwise-averaged transverse power spectra of kinetic energy in mixing layer Li at t = 3.0 ms 
for cases S-32 through S-1. All cases show a turbulent power spectra with a small but significant inertial range, (b) Relative 
power spectra of kinetic energy as a fraction of the ensemble mean for all six cases. While low k power is highly variable and 
shows no clear dependence on Wc, the high wave number portion of the spectra shows a clear dependence on Njc- 


form of the kinetic energy field in the transverse direc¬ 
tions at each layer in the stream wise direction. We then 
sum the 2D spectra along arcs of constant wave number 
magnitude km^ which we define as 

km = ^Jkl^ hi. ( 20 ) 

In this way our transverse spectra for each slice in the 
streamwise direction is designed to avoid preferring ei¬ 
ther the y or z direction. Finally, we average the trans¬ 
verse spectra in the streamwise direction over each mix¬ 
ing layer. Ideally we would like to take a 3D power spec¬ 
tra in order to both study the anisotropy of the transi¬ 
tional turbulence, but for these simulations there is in¬ 
sufficient resolution in the streamwise direction for such 
a procedure. 

Figure [^a) shows the streamwise-averaged transverse 
kinetic energy power spectra for Li at t = 3.0 ms in cases 
S-1 through S-32. All six cases show similar power spec¬ 
tra with energy primarily in the large-scales with ^ 2 
cm“^. Scales with 2 ^ k ^ 6 cm“^ are consistent with 
the establishment of a rough inertial range of turbulence 
and display an approximate scaling as km^^^. The effects 
of numerical dissipation are clearly evident on scales with 
k ^ S cm“^. 

To highlight the variation between cases, we compute 
the ensemble mean spectra by averaging the power spec¬ 
tra over all six cases at each wave number and then di¬ 
viding the power spectra for each case by the ensemble 
mean power spectra. Figure [T^b) shows the fraction of 


kinetic energy power per mode for each case relative to 
the ensemble mean. We have purposefully omitted modes 
with km > ^ cm“^ to focus on the energy-containing and 
inertial scales of the flow. The energy-containing scales 
are quite noisy and do not reveal any clear dependence 
on the initial spectral density of the simulation, however 
the inertial scales do. For 2 < k < 6 cm“^ the power 
per mode is an almost monotonic function of the initial 
spectral density. Cases with lower initial spectral density 
show generally more power in their inertial scales, indi¬ 
cating that power has been less efficiently moved to small 
scales where it can then be dissipated. 


V. SUMMARY AND CONCLUSIONS 

In this paper we have presented a series of numerical 
experiments investigating shock-driven turbulent mixing 
due to the Richtmyer-Meshkov instability. Our models 
include both rarefaction and reshock, which trigger tran¬ 
sition to turbulent behavior. Specifically, we have in¬ 
vestigated the effect of the initial spectral density of the 
perturbations, or the number of Fourier modes over a 
given wave number band, applied to the two material 
interfaces as IC. We have shown that increasing the ini¬ 
tial spectral density leads to enhanced mixing, especially 
after reshock when the mixing layers show more well- 
developed turbulence. We have also shown that evidence 
from our simulations indicates that the primary cause of 
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this variability in mixing is the inhibition or promotion 
of the turbulent cascade. When the initial spectrum of 
perturbations is sparse additional time is required for the 
turbulent cascade to populate the Fourier modes needed 
for the establishment of an inertial range, while a dense 
initial spectrum of interface perturbations is much more 
readily adapted to provide the needed turbulent cascade 
to small scales where motions and material can be dif¬ 
fused by our numerical treatment. 

We have also investigated the dependence of initial 
spectral density on the width of the mixing layers at 
both the front and back of our high-density band. In 
both cases no clear dependence on initial spectral den¬ 
sity is observed. This indicates that the initial spectral 
density has little impact on the formation and growth of 
bubbles and spikes for broadband initial perturbations. 

More broadly, this work reinforces the growing under¬ 
standing of the importance of IC in determining the late¬ 
time mixing properties of RMI-driven turbulence. Previ¬ 
ous theoretical, computational, and experimental efforts 
have explored the impact of single and multi-mode char¬ 
acteristics, spectral slope, and amplitude of initial pertur¬ 
bations [e.g., mini El ESI EH]. To this we presently add 
initial spectral density as yet another parameter which 
shows clear impact on the mixing behavior of RMI-driven 
turbulence. 

Finally, we note that this work may be of particular 
interest to the development of modal models designed 
to initialize RANS turbulence models [49]. Such mod¬ 
els must specify some set of modes to track, however a 
detailed assessment of best practices for specifying this 


initial set of modes remains to be done. Our work demon¬ 
strates the importance of the number and density of ini¬ 
tial modes in 3D and also introduces the shocked heavy 
band problem. We envision the shocked heavy band 
problem could be a relevant test case to assess initializa¬ 
tion procedures and benchmark predictions with current 
reduced-dimension (ID, 2D) RANS of turbulent mixing. 

In conclusion, we have demonstrated that for shock- 
driven turbulent mixing the choice of spectral density of 
the initial perturbations has a clear effect on the mixed¬ 
ness achieved - particularly at late times after reshock. 
Our analysis indicates that this effect is largely the re¬ 
sult of an inhibited turbulent cascade for cases with low 
initial spectral density. In contrast, the total kinetic en¬ 
ergy and mix width were not significantly impacted by 
changing the initial spectral density. 
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